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We describe an example of a structurally stable heteroclinic network for which nearby orbits exhibit 
^irregular but sustained switching between the various sub-cycles in the network. The mechanism for 
^^^witching is the presence of spiralling due to complex eigenvalues in the flow linearised about one of 
Q^the equilibria common to all cycles in the network. We construct and use return maps to investigate 
^^jthe asymptotic stability of the network, and show that switching is ubiquitous near the network. 
ly^Dome of the unstable manifolds involved in the network are two-dimensional; we develop a technique 
•to account for all trajectories on those manifolds. A simple numerical example illustrates the rich 
' dynamics that can result from the interplay between the various cycles in the network. 

o> : 
o ■ 

J>l| Introduction 

^^Heteroclinic cycles and networks are invariant sets that can occur in a struc- 
■ — turally stable way in systems with symmetry, and are known to provide a 
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robust mechanism for intermittent behaviour in these systems. For the pur- 
poses of this paper, we adopt the following definitions of heteroclinic cycles 
and heteroclinic networks. In the literature, there are more complicated defi- 
nitions [6], but the simpler definitions presented here suffice for our purposes. 
For a finite-dimensional system of ordinary differential equations (ODEs), we 
define: 

Definition. A heteroclinic cycle C is a finite collection of equilibria 
{£i,---,£n} of the ODEs, together with a set of heteroclinic connections 
{71 (t), . . . ,7„(t)}, where Jj(t) is a solution of the ODEs such that jj(t) — > £j 
as t — > —00 and jj(t) — > as t — > oo, and where Cn+i = £1. 

Definition. Let Ci,C2, ... be a collection of two or more heteroclinic cy- 
cles. We say that M = (Ji^-i forms a heteroclinic network if for each pair 
of equilibria in the network, there is a sequence of heteroclinic connections 
joining the equilibria. That is, for any £j,£k G M, we can find a sequence of 
heteroclinic connections {7 Pl (t) , . . . , 7 P , (t) } G J\f and a sequence of equilibria 
{Ur • • >£m !+1 } G such that £ mi = £j, £ mi+I = and 7 Pi is a heteroclinic 
connection between £ m . and ^ m<+1 . 

Under this definition, a heteroclinic network is a connected collection of 
heteroclinic cycles, possibly infinite in number. We allow for an infinite number 
of cycles to co-exist in a network, as can occur, for instance, when one of the 
equilibria has a two-dimensional unstable manifold and there is a continuum 
of heteroclinic connections between that equilibrium and another. However, 
we restrict to the case where the set of all equilibria in the network is finite. 
In general, heteroclinic orbits can connect invariant sets other than equilibria, 
such as periodic orbits or chaotic saddles; we do not consider this possibility 
here. 

Structurally stable heteroclinic cycles in symmetric systems have been stud- 
ied extensively in recent years, with a canonical example arising in the context 
of rotating Rayleigh-Benard convection [14] being analyzed in [19]. A good 
deal is now known about conditions for existence and stability of heteroclinic 
cycles (e.g., [9,22,24,26]), and some results are also known about bifurcations 
of heteroclinic cycles and networks (e.g., [10,16,17,31,35]) and the effect on the 
dynamics of small symmetry breaking (e.g., [15,20,25,27,34]) or the addition 
of noise (e.g., [3,12,36]). Some experimental observations of near-heteroclinic 
cycles have been reported (e.g., see [28] for a recent example); experimental 
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noise and small symmetry-breaking effects prevent exact heteroclinic cycles 
from occurring, but near-heteroclinic structures are seen in certain regimes. 

The dynamics near networks of heteroclinic cycles has been studied in, for 
instance, [1,5,13,17,21,30]. There are some natural questions to ask about the 
dynamics near heteroclinic networks. For example, is it the case that one cy- 
cle in the network is dominant, in that most trajectories near the network are 
attracted to that cycle, resulting in the network structure not being observed? 
Can more than one cycle be observed? Are there trajectories that switch be- 
tween the cycles in the network in a sustained way, visiting all parts of the 
network eventually? Partial answers to these questions have been established. 

Krupa and Melbourne [23] find conditions under which one cycle dominates 
in certain cases. However, their analysis does not cover the example of interest 
in this paper. 

Kirk and Silber [21] construct an example where more than one cycle is 
observable, specifically showing that open sets of trajectories near the network 
may be attracted to each of the two primary cycles in their network. Despite 
both cycles being observable, there is no sustained switching in this example: 
an orbit may switch from one cycle to the other initially but may not switch 
back again. The effect of small noise on the network in [21] is studied in [3], 
where it is shown that noise can either induce switching of trajectories between 
cycles in the network or enhance the attractivity of certain cycles. 

Aguiar et al. [1] describe an example, motivated by conjectures of Field [18] , 
where trajectories switch between excursions about different cycles in a het- 
eroclinic network. In this example some of the connections in the network 
result from transversal intersections between stable and unstable manifolds of 
equilibria, with the consequence that the network probably does not attract 
open sets of initial conditions. Because of the reinjection mechanisms built into 
the network, trajectories will make repeated passes near the transversal inter- 
sections and there are trajectories that follow arbitrarily complicated paths 
around the network, but these trajectories will not approach the network. 
Another example of this type is found in [2] . 

Postlethwaite and Dawes [30] examined an example of a heteroclinic network 
in which trajectories can exhibit periodic or aperiodic patterns of excursions 
past the various cycles in the network. The mechanism inducing switching 
between cycles in this network is a transverse instability of each cycle in one 
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direction. Almost all trajectories near a cycle eventually leave that cycle for 
another cycle, but that cycle in turn is unstable in a transverse direction and 
trajectories eventually leave that cycle too, ultimately returning to a neigh- 
bourhood of the original cycle. This mechanism operates when the network 
as a whole is essentially asymptotically stable, so trajectories get closer to the 
network as they cycle around the network. 

Ashwin et al. [12] describe irregular switching near a heteroclinic network 
connecting periodic orbits and chaotic saddles confined to two invariant sub- 
spaces. In this case, the mechanism determining the switching is said to be 
nonlinear, since it appears to operate in a part of phase space well away from 
the invariant subspaces containing the periodic orbits and chaotic saddles. 
In this example, structurally stable connections again arise from transversal 
intersections of manifolds. 

In this paper, we present another example in which orbits near a hetero- 
clinic network switch repeatedly between excursions about the different cycles 
in the network. The mechanism for switching is the presence of a pair of com- 
plex eigenvalues in the linearisation of the flow about one of the equilibrium 
solutions in the network. Unlike the examples in Aguiar et al. [1,2], where 
complex eigenvalues also occur, the heteroclinic connections in our network 
are structurally stable because of the symmetries of the problem and are non- 
transversal, with the consequence that, so long as symmetries are preserved, 
the network can attract open sets of initial conditions. As discussed further 
below, our example can exhibit an interesting form of switching, where the 
network structure is evident in the long term dynamics even though the net- 
work is not attracting. There are similarities between this and the switching 
observed in [1,2], as discussed below. 

The network we study is in R 4 with Z| symmetry, and is shown schemati- 
cally in Figure [TJ It consists of a set of six equilibria denoted A, B, X, Y, P 
and Q, their conjugate copies under action of the symmetries, and the set of 
heteroclinic connections joining the equilibria. Some of the heteroclinic con- 
nections occur in two-dimensional families, as indicated in Figure [TJ There are 
many different heteroclinic cycles evident in Figured! e.g., A — > B — > X — > A 
and A — > B ^ P ^ X ^ A. The network is the union of all these cycles, and 
is described in more detail in section [2j 

Our analysis of this example proceeds in a standard way via construction 
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Figure 1. Schematic diagram showing part of the heteroclinic network studied. For clarity, the 
equilibrium Q is not shown; this equilibrium plays a similar role to equilibrium P except that the 
one-dimensional heteroclinic connections from Q connect to —X and Y instead of X and Y. The 
remaining (conjugate) parts of the network are obtained under the action of the Z| symmetry 
group. The thin curves represent single heteroclinic connections while the bold curves indicate that 

a two-dimensional family of connections exists between the relevant equilibria. The double 
arrowhead on the connection from B to X indicates that expansion near B in the direction of this 
connection is stronger than the expansion in the direction of the connection from B to Y. 

of return maps that approximate the dynamics in a neighbourhood of the 
heteroclinic network. A feature that complicates the construction is the exis- 
tence of two-dimensional unstable manifolds of some of the equilibria in our 
network and hence of continua of heteroclinic connections between some pairs 
of equilibria. A novel aspect of our work is the way in which we allow for 
this complication; we have developed a relatively straightforward way to keep 
track of trajectories near two-dimensional unstable manifolds even when differ- 
ent orbits within a manifold connect different pairs of equilibria. The method 
is related to the technique used in [33] to analyze a homoclinic bifurcation 
in a problem with a two-dimensional unstable manifold, and it extends pre- 
vious work on other problems with two-dimensional unstable manifolds such 
as [4,5,8,9,20]. 

Using these techniques, we are able to find a simple condition for asymptotic 
stability (resp., instability) of the network. The condition is as expected: the 
product of the ratio of contracting to expanding eigenvalues seen by a tra- 
jectory as it traverses the network must be greater than one (resp., less than 
one), regardless of the itinerary of the trajectory though the network. There is 
also an intermediate case, where whether there is net contraction or expansion 
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depends on the itinerary of the orbit past the various equilibria in the network. 
We go on to show that switching is ubiquitous in the network, whether or not 
the network is asymptotically stable. If the network is asymptotically stable, 
then we show that while orbits generically continue to switch as they approach 
the network, visits to certain equilibria become increasingly rare. 

We find that a particularly interesting form of switching can occur in our 
network. If one of the cycles within the network attracts trajectories (i.e., a 
trajectory ends up closer to the network after making one passage near that 
cycle) while other cycles repel trajectories, then the net effect can be that a 
typical trajectory approaches an attractor (possibly chaotic) that lies near the 
network, with the trajectory repeatedly (but not uniformly) passing close to 
all parts of the network even though the network is not itself attracting. Under 
this scenario the network structure will be observed in the long term dynamics 
even though the network is not attracting. We report numerical observations 
of this form of switching, and defer a detailed investigation to a future paper. 

The rest of this paper is organised as follows. $2] contains a description 
of our heteroclinic network and details of construction of the maps used to 
approximate the dynamics near the network. In $3] we find a condition for 
asymptotic stability of the network, and derive results about switching near 
the network. £0] gives results from numerical simulations of a system of four 
ordinary differential equations, illustrating the various switching phenomena 
associated with our example. Conclusions are contained in $5) 



2 The heteroclinic network 

We consider a system of ordinary differential equations, x = f(x), where 
x = (xi, X2, X3,ys) £ K 4 and f : R 4 — > R 4 is a C 1 vector-valued function. We 
assume this system is Z^-equivariant with the following equivariance proper- 
ties: 




i = 1,2,3 



(1) 
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Ki : {x 1 ,x 2 ,x 3 ,y 3 ) -> (-x 1 ,x 2 ,x 3 ,y 3 ), 
k 2 : {x 1 ,x 2 ,x 3 ,y 3 ) -> (xi,-x 2 ,x 3 ,y 3 ), 
k 3 : (xi,x 2 ,x 3 ,y 3 ) -> (x 1 ,x 2 ,-x 3 ,-y 3 ) 



These symmetries ensure the existence of dynamically invariant subspaces in 
which robust saddle-sink connections can occur. We make the following as- 
sumptions about the dynamics in these subspaces (see Figure [2]). 

• Al: There exist symmetry-related pairs of equilibria ±A and ±B on the x\ 
and x 2 coordinate axes, respectively. Within the invariant plane x 3 = y 3 = 0, 
A is a saddle and B is a sink and there is a heteroclinic connection from A 
to B. See Figure da). 

• A2: There exist symmetry-related pairs of equilibria ±X, ±Y, ±P and ±Q 
in the invariant plane x\ = x 2 = 0. Within this subspace, ±X and ±Y are 
sinks, while ±P and ±Q are saddles. The eight equilibria together with the 
heteroclinic connections between them make up an invariant curve C, which 
is topologically a circle. We hereafter refer to C as a circle, and we assume 
that C can be parametrised by the angle 6 3 , the polar angle in the {x 3 ,y 3 )- 
plane. Note that the intersections of the stable manifolds of ±P and ±Q with 
the invariant plane form the boundaries between the basins of attraction of 
±X and ±y in the invariant plane. Only a small part of each intersection 
is shown in Figure E^b), to avoid giving a misleading impression about the 
dynamics near the origin of the (x 3 , y 3 )-plane, but each intersection curve 
in fact extends to the origin of the subspace. 

• A3: Within the invariant subspace x\ = 0, there exist two-dimensional 
manifolds of saddle-sink connections from B to ±X and ±Y (Figure [2(c)). 
There are also one-dimensional (saddle-saddle or saddle-sink) heteroclinic 
connections from B to dbP and ±Q and from ±P and ±Q to ±X and ±Y, 
as shown in Figure E^c). The unstable manifold of B is two-dimensional, 
and the stable manifolds of ±X and ±Y are each three-dimensional within 
the subspace. 

• A4: Within the invariant subspace x 2 = 0, there exists a two-dimensional 
manifold of saddle-sink connections from ±X, ±Y, ±P and ±Q to A. 
Within this manifold, A is a stable focus. A similar manifold connects the 
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(c) (d) 

Figure 2. Dynamics within the subspaces invariant under the symmetries K\, K2, ft3 and their 
combinations. For clarity, only part of the relevant subspaces are shown in panels (a), (c) and (d), 

with the dynamics in the omitted parts being obtained by applying the symmetries, (a) The 
invariant plane X3 = yz = 0, showing the heteroclinic connection from A to B. (b) The invariant 
plane x\ = 1E2 = 0, showing the invariant circle C and the equilibria ±X, ±V, ±P and ±Q that lie 
on C. (c) The subspace xi = showing part of the two-dimensional unstable manifold of B and 
part of the circle C in the (2:3, 1/3) plane. Most trajectories leaving B go to either ±X or iY, but 
some isolated trajectories go to ±P or ±Q. For convenience, the equilibria ±X and ±Y are chosen 
to lie on the coordinate axes, with the eigenvectors of the corresponding linearised flow at B aligned 

with the axes, but they are not constrained by symmetry to be there, (d) The subspace X2 = 
showing spiralling of the unstable manifolds of X , Y and P into A. The unstable manifold of Q (not 
shown) behaves similarly. In each subspace, the flow is strongly contracting in the radial direction. 

equilibria on C to —A. Apart from the heteroclinic connections from ±P 
and ±Q to ±X and ±V, the unstable manifolds of ±P and ±Q are con- 
tained in the stable manifolds of A and —A. There are no equilibria other 
than the origin and those mentioned above lying in the subspace X2 = 0. 
See Figure H^d). 
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• A5: Equilibrium B has real eigenvalues corresponding to dynamics in its 
unstable manifold, and these eigenvalues are unequal. 

Assumptions Al— A5 ensure the existence of the heteroclinic network shown 
in Figure [H The symmetries k% and ki ensure that x\ and X2 cannot change 
sign along a trajectory, so we consider x\ > and X2 > only. The complex 
eigenvalues at A enable both signs of X3 and 2/3 to occur along trajectories. To 
simplify our analysis, we make the further assumptions: 

• A6: At A and B, the directions of strongest contraction lie along the coor- 
dinate axes x\ and X2 respectively. At each of ±A, ±Y, ±P and =LQ, the 
direction of strongest contraction lies in the (£3,2/3) plane; this direction is 
automatically transverse to the connections from ±P or ±Q to ±X or ±Y. 

• A7: The two expanding eigenvectors at B lie in the X3 and 2/3 directions. 
Without loss of generality we assume that the eigenvalue in the X3 direction 
is larger than that in the 1/3 direction. We also assume that the linearisation 
around A, where there are complex eigenvalues, is in Jordan form. 

• A8: The equilibria ±X and ±Y are, respectively, on the X3 and y^ coordi- 
nate axes. 

Note that we can always choose coordinates so that at least one of A7 
and A8 is satisfied, but we assume both are satisfied in order to simplify the 
calculations. This has no effect on our results. 

Thus the overall network is A — > B — > C — > A, where, within C, trajectories 
can visit any of ±A, ±Y, ±P and ±Q, although only in certain orders as 
indicated in Figures Q] and [2 All cycles in the network contain either three or 
four equilibria. 

Definition. For a trajectory <p(t) close to the network, we define the 
itinerary of <fi(t) to be the sequence {^} of equilibria visited. That is, 
is the itinerary of the trajectory <j){t) if there exists an increasing sequence of 
times {tj} such that the distance from <f>(tj) to the equilibrium £j is less than 
some small constant. For a trajectory that stays close to a single heteroclinic 
cycle, the itinerary will be a periodic sequence, with the (minimal) length of 
the repeating segment of the sequence being equal to the number of equilibria 
in the cycle. 
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In this paper we use the following definition for switching, defined for a 
particular trajectory close to a network. Note that Aguiar et al. [1] define 
switching as a property of a network, not individual trajectories. 

Definition. We say a trajectory switches if, as t — > oo, the itinerary is not 
eventually a periodic sequence with minimal period three or four, that is, the 
trajectory does not eventually remain near a single cycle of the network. In 
this definition, we distinguish between conjugate equilibria, that is, X and — X 
and so on. 

2.1 Coordinates, cross-sections, and local maps 

In this section, we define the coordinates, cross-sections, and local maps re- 
quired for modelling the dynamics near our heteroclinic network. 

Near A and B, we define local coordinates that place the equilibrium at the 
origin. Assumption A7 guarantees that the coordinate axes are aligned with 
the eigenvectors of the relevant linearised system. We use polar coordinates 
when it is more convenient: (0:3,2/3) becomes {r 3 ,9 3 ), where X3 = r3Cos#3 and 
2/3 = r 3sin9 3 . We write Xi or if the local coordinate is the same as the 
corresponding global coordinate. At A and B we use u\ and u 2 for the radial 
coordinate relative to the equilibrium points. At the invariant circle C, we use 
a ^-dependent coordinate transformation to define a radial coordinate U3. 

Near A, the linearised flow is given by: 

ui = -taui, &2 = e A x 2 , x 3 = -c A x 3 - uy 3 , y 3 = ux 3 - c A y 3 , (2) 

where r A , e A , c A and uj are positive constants. In polar coordinates, the x 3 
and y 3 equations give r 3 = —c A r 3 and 6 3 = lo. 
Cross-sections near A are defined as: 

= {(u 1 ,x 2 ,r 3 ,9 3 ) I \ux\ < h,0 < x 2 < h,r 3 = h,0 < 3 < 2tt}, 
H^ ut = {(«!, x 2 , r 3) 9 3 ) \\ Ul \<h,x 2 = h,0<r 3 <h,0<9 3 < 2ir}. 

Here < h 1 is a parameter small enough that the cross-sections lie within 
the region of approximate linear flow near A (and similarly near B and C, as 
required below). 

The flow near A induces a map 4> A : Hjj -> H^ ut , which is obtained to 
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lowest order by integrating equations (|2|): 



Mui,x 2 , h, e 3 ) = L (^) 7 \h,h (^) 5a ,e 3 - ^ log (^) \ (4) 



CA 
eA ' 

Near B, the linearised flow is: 



where 8 a = —■ 



XI = ~C B Xi, U 2 = -TBU2, x 3 = e Bx X3, 2/3 = eByV3, (5) 

where rg, esx-, &By, cb ar e positive constants. Prom A7, we have eg x > &By 
Cross-sections near B are defined as: 

= {(x 1 ,u 2 ,r 3 ,9 3 ) \xi = h, \u 2 \ < h,0 < r 3 < h,0 < 6 3 < 2vr}, 
H£ ut = {(xx,u 2 , r 3 , 63) I < Xl < h, \u 2 \ <h,r 3 = h,0<8 3 < 2vr}. 

The flow induces a map <f>s '■ H 1 ^ — > H B nt , which is obtained to lowest order 
by integrating equations (J5]). The map cannot be written down explicitly, but 
is computed as follows. First, the x 3 and y 3 equations are solved: 

x 3 (t) = r 3 (O)cos0 3 (O)e e ^, y 3 (t) = r 3 (0) sin 6 3 (0) e eB -\ 

where r 3 (0) and 3 (O) are the initial values of the radial coordinates (i.e., on 
H"^). The trajectory crosses H^ 1 * when r 3 (t) = h, so the transit time Tb is 
found by solving the equation 



h \ 2 

= cos 2 3 (0) e 2eB * TB + sin 2 3 (O) e 2eB ^ B (7) 



for Tb in terms of r 3 (0) and 3 (O). This yields the local map 4>b '■ — > H^ 1 *: 

4> B {h,u 2 ,r 3 ,e 3 ) = (he~ CBTB ,u 2 e~ rBTB ,h,tan^ (tan(0 3 )e (e ^- e ^ )Ts )) . (8) 



For later convenience, we define S Bx = <5 mm — and 8b v = 5"^ ax = 

The treatment of the dynamics near the invariant circle C is more compli- 
cated. We assumed in A2 that C can be parameterised by the angle 9 3 . The 
rate of relaxation onto C is controlled by the 3 -dependent eigenvalue — rc{0 3 ). 
The assumption of strong contraction in the radial (r 3 ) direction (A6) means 
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that the dynamics on C of 9 3 can be described by a one-dimensional nonlinear 
ODE of the form #3 = g{9 3 ). The presence of ±X and ±V on the coordinate 
axes will require g(0) = g(ir/2) = g{it) = g(37r/2) = 0. The presence of ±P 
and ±Q will require further zeroes of g. This results in the flow near C being 
given by: 

xi = &c{Q?,)x\, x 2 = -cc{0 3 )x 2 , u 3 = -rc(03)u3, O3 = 3(6*3), (9) 

where rc, ec and cq are positive functions of 6*3. 
Cross-sections near C are defined as: 

= {(xi,x 2 ,u 3 ,9 3 ) |0 <xi < h,x 2 = h, \u 3 \ < h,0 < 6 3 < 2tt}, 
H£ ut = {(xi, x 2 , u 3 , 9 3 ) I x x = h, < x 2 < h, \u 3 \ < h, < 63 < 2vr}. 

There is a continuum of heteroclinic connections from B to the various equi- 
librium points in C, and defining the cross-sections in this way allows us to 
keep track of all these connections. 

The flow induces a map (fic '■ — > H^*. As in the case of the flow past B, 
we cannot write down the map explicitly, but it is computed as follows. First, 
the 6*3 equation is solved using an initial condition #3(0), yielding 6*3 (t). Then 
the x\ and x 2 equations are solved: 

zi(t) = xi(0)exp^ e c {93{t'))dt^j, x 2 {t) = hexp (- cc(6 3 (t'))dt' 

The trajectory crosses when xi(t) = h, so the transit time Tc can be 

found in principle by solving 

\ C (9 3 (t'))dt> = -log (^) (11) 

for Tq in terms of the initial values xi(0) and 6*3(0) on H™. Then the local 
map 4>c '■ H™ — > H^* is given by 

<t>c(xi,h,u 3 ,e 3 )= (h,hexp(-^ C c c (03(t'))dt'^ , u 3 (T c ) , 9 3 (T C )^ , (12) 
where U3{Tq) = 1*3 exp (- J Tc r c (0 3 (f)) dt') . For later convenience, we define 
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$cx, 5cy, $cp and 5qq, to be the ratio ^(63) evaluated at the points X, Y, 
P and Q, respectively. 

Neither of the local maps (f) B and (j>c can be written down explicitly. In 
the case of (j) B , the obstruction is only that we cannot write down an explicit 
solution of (J7J) for the transit time T B . In the case of 4>c the nonlinear evolution 
of #3 within C cannot be written down explicitly. However, in both cases, 
we will be able to give bounds on some properties of the trajectories, and 
this turns out to be sufficient for the purposes of determining stability and 
switching properties of the network. 

2.2 Global maps 

We construct global maps Vl/jj to approximate the dynamics near the hetero- 
clinic connections between A, B and C. In each case, we linearise the dynamics 
about the unstable manifold leaving the invariant set, taking into account the 
fact that the unstable manifold of A is one-dimensional, but the unstable man- 
ifolds of B and C are two-dimensional. We make use of the equivariance of the 
vector field in our map construction. 

The simplest of the global maps is ^ ab ■ H^ 1 * — > H"J. The heteroclinic 
connection from A to B intersects H^ ut at (u±, x^, £3, f/3) = (0, h, 0,0), and 
intersects at (x±, U2, x 3 , 1/3) = (h, ee,0, 0), for a small constant e#. Gener- 
ically, €b 7^ and we assume that this is the case in the following. Near the 
heteroclinic connection the map expressed in cartesian coordinates is at low- 
est order an affine linear transformation. Converting to polar coordinates, this 
yields, at leading order: 

V AB (u u h,r 3 ,e 3 ) = (h,e B ,D B (9 3 )r 3 ,9 B (6 3 )), (13) 

where D B (0 3 ) is an order-one function of #3 that indicates how the small 
variable r 3 is scaled in the transition from A to B, and 9 B {9 3 ) is an order-one 
function of #3. The invariance of this map under the symmetry k 3 ensures that 
there is no constant term in the r 3 component. The overall effect of this map is 
to multiply r 3 by an order-one amount D B , and to rigidly rotate the angle #3. 

The unstable manifold of B is two-dimensional; it intersects H B nt at 
(xi,U2,r 3 ,9 3 ) = (0, 0, h, 9 3 ), for < 9 3 < 2tt, and it intersects Hg! at 
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(x±, x 2 , Us, 63) = (0, h, €0(63), 9c(6z))i where ec is a small function of 83 and 
9c is an order-one function of 83. For small x\ and u 2 , we have at leading 
order: 

*Bc(xi,VQ,h,6 3 ) = {Dc^xuhec^dcWs)) , (14) 

where Dc(03) is an order-one function of #3. Here ec($3) plays a similar role 
to the constant e# in fjl3[) . except that it takes on a different constant value 
for each heteroclinic connection and so is a function of #3. As with e#, €0(63) 
is generically non-zero and we assume that ec(#3) 7^ for any #3. 

The effect of (|14|) is to multiply the small variable x\ by an order-one function 
of #3, and to map the outgoing angle #3 to an incoming angle Oq. Unlike in 
the case of \£ab> the effect of 9q need not be a rotation. 

For the global map ^ca ■ H^* — ► H™, we also use ^3,63) rather 
than (#3,2/3). The unstable manifold of C is two-dimensional; it intersects 
Hg! 11 at (xi,x 2 ,u 3 ,9 3 ) = (/i, 0, 0, 6> 3 ), where < 63 < 2ir, and it intersects H™ 
at (in, X2, T3, ^3) = (ca(^3), 0, h, 9a(9s)), where €a is a small function of #3, 
and 8a is an order-one function of 83. For small x 2 and 113, we have: 

*CA(h, x 2 ,u 3 ,e 3 ) = (eA(93),D A (9 3 )x 2 ,h,8 A (83)) , (15) 

where Da{93) is an order-one function of #3. The effect of this map is to 
multiply the small variable x 2 by an order-one function of #3, and to map the 
outgoing angle 6*3 to an incoming angle 9 a- As in the case of *sc, the effect 
of 9a need not be a rotation. 

3 Analysis of the maps 

By composing the local and global maps in an appropriate order, we construct 
return maps that approximate the dynamics near the cycles in our network. 
We are interested in finding conditions under which the network as a whole is 
attracting, and in describing the switching properties of trajectories as they 
travel around the network. We are particularly interested in trajectories that 
repeatedly visit both X (or —X) and Y (or —Y). 
In our analysis, as in [20], we make use of the observation that at each cross- 
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section, the four variables play distinct roles. Two variables are unimportant: 
the one that is equal to h, and the one in the radial direction. Of the other two 
variables, one is small and measures the distance from an invariant subspace, 
and the other is the order-one angle #3. At each cross-section, the roles change, 
but there are always small and angle variables. 



3.1 Stability results 

To show that the network as a whole is attracting, we must find conditions 
under which the small variable decreases each time around the network. In 
order to do this, we bound this variable over all possible values of the angle 
variable. This means that we need to take into account the details of which 
part of C is visited by the trajectory. We are unable to compute the stability 
result by direct computation of a return map, since the local maps <pB and 
(j>C are only known implicitly, but the lengthy computation below achieves the 
same result. 

We start on at (u±, X2, h, 83), and consider the effects of maps (pA, ^AB, 
4>B-, ^bc, 4>c and ^ca in turn. We assume that xi <S 1. After cpA and ^ab, 
we arrive on H 1 ^ at: 

w*. (I)) 

where we have discarded all higher-order corrections. For convenience, we label 
the values of and #3 on as r 3 11 and #3 , with a similar convention 
on other cross-sections. We define D^ ax and Dg m to be the maximum and 
minimum values of Db taken over all values of #3. Then we can bound the 
small variable r^ Bia ^ by 

hDr{^) SA <r!f-><kDr"(^y". (16) 

Next, we consider the effect of maps <j>B and ^bc- Recall from (A7) that 
e-Bx > ZB y . From ([7]) we can bound Tg: 



< T B < 
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As a result, the small variable x^ 8 ™*' on H B ut is bounded by 
Note that 8b x < 5 By After \P bc, trajectories enter at: 

[d c {ef OM) ) 4 B ^\h,e C (e 3 ),e c 

where we have discarded all higher-order corrections. We define D™ ax and 

Dq 1u to be the maximum and minimum values of Dc taken over all values of 

(C- ) 

#3. Then we can bound x\ '" by 

ftC ^— < xp } < ^— - . (17) 

Finally, we consider the effects of maps 4>c and ^ca- The first of these is 
the most complicated as trajectories can enter the neighbourhood of C close 
to any of the equilibrium points ±X, ±Y, ±P or ±Q (or in between), and can 
similarly exit the neighbourhood of C close to any of the equilibrium points 
±X, ±y, ±P or ±Q (or in between). We take all possibilities into account 
and derive a bound on the exit value of the small variable xi- 

To simplify the discussion, we consider in detail only the case of trajectories 
arriving at C between X and P; the other cases can easily be deduced from 
this one. Within C, X is stable, with a stable eigenvalue —Ax, and P is 
unstable, with an unstable eigenvalue Ap, with Ax, Ap > 0. The corresponding 
eigenvectors are within the (#3, 2/3) plane, and are tangent to C at X and P. To 
aid the analysis, we consider cross-sections orthogonal to C, at 63 = 9f +h and 
63 = 9^ — h, where 6^ and 9^ are the #3 coordinates of X and P respectively, 
and h is as before (see Figure [3|). Note that we have defined our coordinates so 
that 63 = 0, but for clarity of the following discussion we leave this constant 
in symbolic form. 

There are three main possibilities, indicated in Figure [3) (a) the trajectory 
crosses H™ near X and must therefore cross H^* near X as well; (b) the tra- 
jectory crosses H™ near P and also crosses H^* near P; and (c) the trajectory 
crosses near P and leaves the neighbourhood of P, heading towards X, 
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(c) (b) 
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Figure 3. Schematic showing three main possibilities of how trajectories pass through the region 
near C. Trajectories are shown as dashed lines, Poincare sections as solid bold lines. Equilibria are 
indicated by dots. The labels (a), (b) and (c) identify trajectories representative of three cases 

discussed in the text. 



and so crosses H^* near X. There are the additional possibilities that the 
trajectory crosses or H^* in between X and P\ we discuss these cases 
below. Throughout this discussion, we disregard the radial coordinate r$. 

Case (a) is straightforward: the flow near X is given by the linearisation 
of © around X: 

ii = ec(6£) Xl , ± 2 = -cc(eg)x 2 , e 3 = -\x(03 - eg), (is) 

(C- ) (C- ) 

and so the linearised map near X takes an incoming point (x\ , h, 0% m ) to 

where 5c x is the ratio of contracting to expanding eigenvalues evaluated at X. 

Similarly, case (b) is treated by linearising © around P, and results in exit 
values 

where c>cp = Cc 1op1 i s the ratio of contracting to expanding eigenvalues eval- 
uated at P. The condition that the trajectory crosses H^S 1 * before it crosses 
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9^ — h amounts to: 



,(Cn) 



.(On) ' 



< h, 



that is, the trajectory must enter C close enough to P that the small variable x\ 
grows to size h before the angular separation #3 — 6^ grows to size h. 

Finally, in case (c), there are three stages: linearised dynamics near P, a 
jump from P to X, and linearised dynamics near X. The time for the first 
stage is found by setting 9s(t) = 9% — h, and then the {x\,X2) coordinates on 
this section are: 











c c (£»f ) \ 






Xp \ 




h 





Then there is a jump from #3 = 6^ — h to #3 = Q$ + h, during which x\ 
and X2 change by a factor of Dx,i and Z?x,2 respectively. The values of Dx,i 
and Dx,2 depend on h. Lastly, there is the linearised dynamics near X, which 
results in an exit value of X2 given by: 



hD 



X,2 



Dx,ix\ 



h 



(19) 



This case only occurs if 



*3 



> 



(20) 



(otherwise we would be in case (b)). This condition allows us to bound the exit 

M C in) /IP 

values of X2 . Note first that 



h 



< 1. Then if 5cp > <fcx> from ([19 



we have that is bounded above by a constant times I x 1 



ditionally, using (|20j) in (jl9|) . we have that x 



.(Co* 



. Ad- 



is bounded below by a 



constant times 



(4 



(Cm) 



Similar considerations in the case that 5c p < Sex 
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give additional constraints which altogether result in: 

/ (Cm) \ max ( <5 c-p,'5 cx ) / /q\ \ min(5 C p,<5cx) 

hD min I j < x (C out ) < ^max / ^1_ \ 



where Z)™ n and D™ 8 * are constants that depend on Dx,i and -Dx,2 and some 
exponents. 

The same analysis can be used in the cases where trajectories enter or leave 
the neighbourhood of C in between the neighbourhoods of P and X, with 
only minor alterations of the values of the constants D^ m and L>^ ax . Then all 
possibilities (a), (b) and (c) can be assembled, as well as including trajectories 

that visit the equilibria Y and Q as well, and the map ^>ca can be applied. 

(A- ) 

All this results in a bound on the value of x 2 at of the form: 

/ (c , ,„)\ <5 " ax / (c m )\ <5 " ,n 



where we interpret <5™ ax as max (c\7Xj &CYi ^CPi $cq)i and similarly <5™ m , and 
the constants L>™ m and _D™ ax are the smallest and largest of all the constants 
in the local and global parts of the maps. 

Recall that the trajectory started on with particular values of X2 and 
63. In computing these bounds, the value of #3 has been lost, but inequalities 
(fT6j) . (fT7|) and (f21"j) provide the smallest and largest possible values of X2 once 
the trajectory returns to H™: 



hD ABc{f) <x^> <hDTic{-f) ■ (22) 

Here 5 max = 5 A S^ x 5% ax , 5 min = 5 A 6f n 5^ in , and all the constants have been 
amalgamated into D^g C and D™"^ c . 

We have thus established a condition for asymptotic stability or instability of 
the network. If 5 mm > 1, then a trajectory starting close enough to the network 
will return closer to the network (with a smaller value of X2) regardless of which 
itinerary it takes and regardless of the values of D™ 1 ^ and D^gc > an< ^ so the 
network is asymptotically stable. If $ max < l 5 then a trajectory starting close 
to the network will return further away from the network (with a larger value 
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of x<i) regardless of which itinerary it takes, and so the network is unstable. 
The values of the constants can be scaled away in both these cases. 

These conditions for asymptotically stability and instability are as expected: 
the product of the ratio of the contracting to expanding eigenvalues should 
be greater or less than one regardless of the itinerary. The more interesting 
and complicated case is when <5 mm < 1 < <5 max , in which case it appears 
that whether there is net contraction or expansion depends on the itinerary. 
We present some numerical results relevant to this case in section [H These 
results suggest that the network may be essentially asymptotically stable or 
unstable, depending on which of the routes through the network is responsible 
for <5 max > 1 and which is responsible for 5 mni < 1. 

3.2 Switching near the network 

In this section, we show that close enough to the network, there are trajectories 
that, over the course of two circuits around the network, visit any combination 
of the equilibrium points within C in any order. This occurs whether or not 
the network is asymptotically stable. We also show that when the network is 
asymptotically stable, most trajectories repeatedly visit both X and —X as 
they approach the network. On the assumption that the complex eigenvalues 
at A mix trajectories effectively, we estimate how often trajectories visit ±Y 
and show that, when the network is asymptotically stable, visits to ±Y become 
rare as trajectories approach the network. Finally, we show that interesting 
switching dynamics might be possible in the case 5 mm < 1 < <5 max , when some 
parts of the network are attracting and other parts are repelling. 

Figure H] shows schematically how trajectories visit different parts of C, 
according to where they cross H 1 ^. The majority of trajectories go to X or 
— X, and there are cusp-shaped regions that visit P then X, P only, etc., on 
their way to A. 

Consider a line segment of initial conditions on Hg, with a fixed value of #3 
and a range of values of r%: < rjj < r% < r| < h (see Figure [5]). We choose 
63 such that the family of trajectories first travels around the network via the 
point X. We will show that the spread of r$ values translates into a spread 
of #3 values once the line segment has been mapped around the network. This 
arises in particular from the complex eigenvalues at A. We can choose the 
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(a) 

P^X P^Y Q^Y Q^-X 




Figure 4. Schematic diagram showing a slice (with constant 112) of the Poincare section H^, using 
(a) Cartesian and (b) logarithmic polar coordinates. Only part of the slice in (a) is shown in (b). 
Each slice is divided into regions according to the equilibrium in C visited by the trajectories in 
that region as they pass from B to A. 

values of rf and rjj such that after one cycle around the network, the line 
segment covers 63 £ [0, 2n]. This means that there are trajectories with initial 
conditions in the initial segment that visit each of the different parts of C on 
the second time around the network. 
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#3 = 6 3 = TT d 3 = 2lT 

Figure 5. Schematic diagram showing a line segment (solid grey line) on (with 
< 7"3 < T3 < r| < h, shown in logarithmic polar coordinates) and its image (dashed grey line), 
after one cycle around the network. The length of the line segment is chosen such that its image 
covers the full range of values of 6*3. Note that in any particular example, the image could be 
considerably more complicated than a straight line. 

In order to show this, we repeat part of the calculation of section 13.11 but 
starting on instead of H™. For any given initial value of r 3 on H"J, using 
inequalities (|17p and (|2ip . the trajectory crosses H™ with a value of x% that 
satisfies: 

hDfl (£) aS " i3 ~ < zpJ (r 3 ) < hD^ (2) W ' , (23) 

where Dgjf = Z)™ in (Z)™ 11 )' 5 ™"' and D™f = D™ ax (Dg 3 *) 5 ^ . On Hj£, the 

(A- ) 

trajectory has a value of #3 close to that with which the unstable manifold 

of X crosses H™, since we chose our initial line segment such that trajectories 

visited X. This is essentially a constant, but the complex eigenvalues mean 

that 9^°^ depends logarithmically on x^ 1 ^ . We want to choose rf and 

(A. \ 

such that even allowing for the range of values of x\ in (|23p . the local 
map (jU) guarantees that the range of values of 6^ A ° ut ^ covers at least [0, 2tt]. 
This requires that 




where the minimum and maximum in this expression are taken over all possible 
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values of #3 on H 1 ^. This last expression can be rewritten as 



. (r%/h) 5m ~\ 2irc A s , fD, 
log t¥tt^ )> — +5 A logl 



max 
CA 



ri/h) J w \ d ca 

Values of rf and 73 for which this is satisfied can clearly be found, regardless of 
the values of the global constants or the values of 5 mm and 5 max , provided that 
w / 0. Disregarding the global map constants (-D™^ x etc.), the inequality is 
satisfied if we choose r§ > exp (27rcA/w(5 max ). By choosing small enough 
(close enough to the network), the length of the initial line segment — r| 
can be made as small as we wish. 

With this choice of 7-3 and r|, the line segment of initial conditions maps 
to at least a full circle on H^ ut and consequently, using the global map (fT3j) . 
also on H 1 ^, since that map rigidly rotates the angle. Therefore, on their next 
circuit around the network, trajectories from within this family could visit any 
of ±X, ±Y, ±P or ±Q. 

The same argument holds with minor changes regardless of the location 
of the initial line segment, so we conclude that close to the network, there 
are trajectories that visit any equilibrium point in C followed by any other 
equilibrium point in C on two consecutive circuits of the network. Since there 
is freedom in choosing the exact location of the line segment, this argument 
implies the same conclusion can be drawn for open sets of initial conditions. 
Thus we have shown that arbitrarily close to the network, there are open sets 
of orbits that switch from any route around the network to any other route. 

The argument above does not require the network to be asymptotically 
stable, and only refers to two consecutive circuits of the network, and so does 
not demonstrate that typical trajectories will continue to switch for ever as 
they evolve. In the remainder of this section we consider long term switching, 
first in the case where the network is attracting (<5 mm > 1) and then in the 
case where 5 mm < 1 < <5 max . 

The argument above implies that of trajectories starting in a typical ball of 
initial conditions close to the network and first visiting X, fewer than half will 
go on to visit X on their second circuit of the network. In the case 5 mm > 1, 
all trajectories starting in a typical ball approach the network, so their values 
of x^"^ get smaller and smaller, and hence the values of 6^"^ for these 
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trajectories get more and more spread out. Therefore, the argument that the 
chance of visiting X on the next cycle around the network is less than half 
continues to hold; in the limit, we expect that the set of trajectories that visit 
only X and never visit —X has measure zero. Thus, typical trajectories should 
visit both zizX roughly equally, and there is always a chance they could visit 
±Y, ±P or ±Q as well. 

We can estimate the probabilities of visiting different parts of the network 
by computing the proportion of trajectories starting on that will visit 
X or Y on their next time around the network. Consider a circle of initial 
conditions on H 1 ^ with r$ = a (a < h) and < #3 < 2ir. Of those in the first 
quadrant, some will visit X, some will visit P and the rest will visit Y (see 
Figure[5|). The boundaries between the different possibilities are of the form 
X3 = Ky$, where a = eBx/&By > 1> and K is a constant that depends on 
which boundary is being considered. To simplify the discussion, we omit the 
details of those trajectories that visit P, and consider only a single boundary 
that separates trajectories that go via X and those that go via Y. 

The intersection of the boundary with the circle = a can be found by 
solving 

a" = xl + y! = y!(l + KV 3 a - 2 ) 

for 2/3 as a function of a, a and K. For small a and for a > 1, the solution is 
approximately 

y 3 = a (l - ^K 2 a 2a ~ 2 ^ , x 3 = Ka a . 

The proportion of trajectories starting on the circle r 3 = a in that visit Y 
is approximately equal to 2x3 /2na = Ka a ~ 1 /ir for small X3. The same pro- 
portion visit —Y, and the remainder are split equally between X and —X. 

In the case that the network is attracting ((5 mm > 1), the value of r% on 
decreases each time around the network, so a — > 0. Trajectories spend increas- 
ingly long periods of time near A, where the eigenvalues are complex, so one 
would expect that the angle 63 becomes essentially a random variable. In this 
case, the chance of visiting Y or —Y goes to zero, and the chances of visiting 
X or —X will both tend to ^. 
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In the case J max < 1 ; where the network is unstable, trajectories leave the 
neighbourhood of the network and no estimates are possible. 

The intermediate case (<5 mm < 1 < <5 max ) offers the interesting possibility 
that trajectories might maintain an average distance from the network, either 
in a periodic or chaotic fashion. In the latter case, one might expect that an 
average, weighted using the probabilities above, of the contraction around one 
part of the network and the expansion around the other part might lead to 
conditions for the existence of a nearby invariant set. Making a weighted aver- 
age in this way only makes sense if trajectories switch irregularly between ±X 
and ±y, so an existence condition could include a requirement for switching. 
In the case of a periodic orbit, the weighting would depend on the itinerary of 
the orbit. This weighted average would also depend on a, which we interpret 
as the average distance from the network. We give numerical examples of this 
phenomenon below, and defer a detailed analysis of this case to a later paper. 



4 Numerical example 

In this section we present some numerical results based on the following equa- 
tions: 

x\ = x\(l — x\ — Ex^), 

x 2 = x 2 {l-xl-Fxl-Gyl), 

x 3 = 23(1 - Ex\ + Hx\ -x\- Dyf) - uy z x\, 

h = 2/3(1 - Ex\ - Dx\ - yf} + ux 3 x\, 

where D, E, F, G and H are parameters that we vary in our numerics. The 
parameter H controls the relative values of the two expanding eigenvalues at 
the point B. Throughout we assume that H > 0, and that D, E, F, G, (F + 
G)/{D + 1) € (1,3), so that assumptions A1-A8 from section [2] are satisfied. 
The values of all parameters and eigenvalue ratios are given in Table [TJ 

To simplify the presentation, we always choose the parameter D so that 
5c p and 5cq are intermediate between 8c x and 8cy- With this constraint, 
the additional combinations of eigenvalue ratios are 5™ ax = max (5cx, $cy), 
8™ in = mm(5 C x,ScY)- Recall that 

cmax x rmin r rmax r srmaxrmax rmin r rmin rmin 

Ob = °By, 8 B = Obx, 8 = 8a8 b o c , 8 = 6a8 b 8 c . 
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Equilibrium point 


Eigenvalues 


Eigenvalue ratios 


A: (1,0,0,0) 


rA = 2, e A = l, c A = E -1 


5 A = E - 1 


B: (0,1,0,0) 


r B = 2, c B = E - 1, 
e&r = 1 + H, e By = 1 


r _B— 1 

= E — l 


A: (0,0, 1,0) 


r c(^3 ) = 2 > A x = - 1, 

ec(^) = l, cc(6X) = F-l 


X TP 1 


Y- fO "H 
2 . ^u, u, u, 


r^(f) Y ) — 2 Av — D — 1 
e c (9j) = l, cc(eX) = G-l 


<5cy = G 1 




r c (e£) = 2, A P = 2g T i, 
ec(^ P ) = 1, ccr(flf) = " 1 






rcr(^) = 2, A Q = 2g=i, 
ec(^ Q ) = 1, cc(^ Q ) = ft¥ " 1 





Tabic 1. Eigenvalues associated with equilibria of equations (24}, and values of eigenvalue ratios. The 
actual eigenvalues at A (for example) are — r a (radial), — ca =b iiu (contracting) and ca (expanding). The 
eigenvalues -Ax, Ap, —\y and Aq all have corresponding eigenvectors tangent to the invariant circle C 



Parameter 


I 


II 


III 


F 


1.63125 


1.63125 


1.61875 


G 


1.671386719 


1.549316406 


1.671386719 


8 Bx = $f n 


1.25 


1.25 


1.25 


5 By = 5r x 


1.28 


1.28 


1.28 




0.63125 


0.63125 


0.61875 




0.671386719 


0.549316406 


0.671386719 


&CP = &CQ 


0.63496867 


0.57453782 


0.62878055 


Sx 


1.01 


1.01 


0.99 


Sy 


1.10 


0.90 


1.10 


^max 


1.10 


1.0342400 


1.10 


^min 


1.01 


0.87890620 


0.99 



Table 2. 

arc u) — '. 
H)/(E- 



Values of the parameters and the eigenvalue ratios for the three examples. The values in common 
1, H = 0.024, D = 1.02 and E = 2.28. The parameters F and G are set using F = 1 + S x (l + 
■ l) 2 and G = 1 + S Y /(E - l) 2 . 
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We make the further definitions 5x = ^a^Bx^cx and Sy = ^A^By^cV- these 
are the eigenvalue ratios relevant to trajectories that leave B along the x% 
(resp., ys) direction and visit ±X (resp., ±Y). 

To illustrate typical dynamics near the network, we present three examples. 
In these examples, we choose values of 5 close to 1 so that trajectories do not 
approach or leave the network too quickly. The first example has 5 mm > 1, and 
the other two have <5 min < 1 < <5 max , with different choices as to whether it is 
5x or Sy that is less than 1. The specific choices of parameters and eigenvalue 
ratios for the different examples are given in Table [2j 

• Example I: 5 m[n = 5x = 1-01 and 5 max = 5y = 1.10. In this example, 
trajectories approach the network, predominantly switching between X and 
-X. 

• Example II: S mm < 5y = 0.90 and J max > 5x = 1-01. In this case trajectories 
behave much as in Example I, even though the network is not asymptotically 
stable, since almost all trajectories visit ±X most of the time and so are 
most heavily influenced by the value of 5x ■ 

• Example III: 5 min = 5 X = 0.99 and 5 max = 5y = 1.10. Here trajectories 
leave a neighbourhood of the network and end up displaying periodic or 
chaotic switching. 

Care needs to be taken with numerical integration of systems with hete- 
roclinic cycles and networks, because of the potential for rounding errors to 
cause qualitatively incorrect results. We first integrated equations (|24p numer- 
ically using the Bulirsch-Stoer adaptive integrator [32], with a tolerance for 
the relative error set to 10~ 12 for each step. We also rewrote equations (I24p 
using logarithmic variables (log x±, log x^, log r$, 63) instead of (x±, X2, £3, 2/3), 
and integrated the converted equations with a tolerance of 10~ 10 , rising to 
10 -8 for trajectories very close to the network. This enabled us to examine 
whether the numerical methods handle the very large dynamic range of the 
variables without being unduly affected by rounding errors. The two meth- 
ods of computing solutions agree to within the specified tolerance when we 
compute periodic trajectories (in calculations II and III), and they agree to 
within the specified tolerance for times up to about 3000 when we compute 
trajectories very close to the network. Beyond this time, trajectories computed 
by the two methods diverge, but the qualitative behaviour of the trajectories 
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Figure 6. Example I: Phase portrait showing the (2:3,1/3) projection of a single trajectory. The 

trajectory leaves B (at the origin in this projection), goes to C (approximately the circle 
x| + y| = 1) in a more-or-less straight line, travels around C towards ±X or ±Y, then spirals in 
to A (also at the origin) before returning to B. 

is the same. We have confirmed that the results are not sensitive to the exact 
value of the the relative error tolerance that we chose. The results shown in 
the figures below are all computed using logarithmic variables. 

Poincare sections were computed using algorithms from [29]. The nodes 
on the network are all simple equilibria lying within coordinate planes, so the 
numerical issues associated with cycling chaos (chaotic dynamics within nodes 
on the network), as discussed for example in [7, 11] do not arise here. 

4.1 Example I: Trajectories approach the stable network 

In this example, 5 mm = 5x = 1-01 and <J max = § Y = 1-10 are both greater than 
one, and the network is asymptotically stable. The phase portrait and time 
series shown in Figures [6] and [7] correspond to a trajectory started from the 
initial condition xi(0) = 0.01 = h, X2(0) = 1 and £3(0) = 2/3(0) = 10~ 5 , and 
illustrate the occurrence of repeated switching in the transient dynamics. As 
expected from section 13.21 a typical trajectory lying near the network makes 
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Figure 7. Example I: Time series for the trajectory shown in Figure [6] The x\ and X2 plots show 
that the trajectory is approaching the network. The X3 plot shows the repeated switching 
between X and —X, and the 3/3 plot shows that visits to ±Y become increasingly rare as the 
trajectory gets closer to the network. 



repeated switches between X and —X and visits ±Y occasionally, but the 
visits to ±y become increasingly rare as the trajectory gets closer to the 
network. 

Figure [8] shows where the trajectory intersects the Poincare section H 1 ^, 
defined here as x\ = h = 0.01, X2 ~ 1 and < h. After leaving B, the 
trajectory visits ±X or ±Y: the symbols indicate which of the four possibilities 
occurs immediately after the intersection marked. The boundaries separating 
regions of H 1 ^ from which trajectories leave for ±X and ±Y can be clearly 
seen, and are consistent with the results sketched in Figured As the trajectory 
approaches the network, travelling from top to bottom in Figure[8l visits to ±Y" 
are not seen for < 10 -20 or so, although they are in principle possible for 
arbitrarily small r%. Visits to X and —X do not occur in a periodic fashion. 
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Figure 8. Example I: Projection of a Poincare section for the trajectory shown in Figure [6] Black 
astcrixes (blue plusses) indicate that the trajectory visits X (—X) immediately after leaving the 
Poincare section; red crosses (cyan boxes) indicate that the trajectory next visits Y (— Y). The 
single point in a box at the top is the first point in the trajectory, and the ten boxed points at the 
bottom are the final points. This figure illustrates how the trajectory approaches the network; the 
trajectory initially visits ±X and ±Y, but as it approaches the network, visits to ±Y become rare 

while switching between X and —X is persistent. The boundaries of the cuspoidal regions are 
indicative of the boundaries between trajectories that have different routes on their next circuit of 
the network (c.f., Figure [2J. They have been chosen to match the available data for larger T3. 

4.2 Example II: Trajectories approach the unstable network 

In this example, <5 mm < 6y = 0.90 and <5 max > 5x = 1-01, and the network 
is asymptotically unstable. However, trajectories that start close enough to 
the network can still approach the network. For instance, the initial condition 
Xl (0) = 0.01 = h, x 2 (0) = 1, a?3(0) = ito(0) = 10" 40 yields the trajectory 
shown in Figure [9] (below the horizontal line). This trajectory never visits ±Y 
but does switch repeatedly between X and —X while getting closer to the 
network. In contrast, the initial condition xi(0) = 0.01 = h, £2(0) = 1, £3(0) = 
2/3(0) = 10 -20 yields the periodic orbit close to the network shown above the 
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Figure 9. Example II: Poincare section. The trajectory with initial condition x±(0) = 0.01 = h, 
2:2(0) = 1, 2:3(0) = 2/3(0) = 10 — 40 (below the horizontal line) approaches the network, starting at 

the single boxed point just below the line, and ending at the ten boxed points at the bottom of the 
figure. The trajectory with initial condition 2:1(0) = 0.01 = h, 2:2(0) = 1, 2:3(0) = 2/3(0) = 10~ 20 

approaches a stable periodic orbit represented by the two boxed points above the line. Symbols are 

as in Figure [8] 

horizontal line in Figure [H There are other stable periodic orbits further away 
from the network. 

This behaviour is consistent with the discussion in section [3j Since 8x > 1, 
a trajectory that starts close enough to the network will mostly only visit ±X 
(almost never ±Y) and so can approach the network even if 5y < 1- On the 
other hand, the discussion in section [3] predicts that there are trajectories 
arbitrarily close to the network that visit ±Y sufficiently often to be repelled 
from the network, and so the network is unstable. We have found examples of 
stable periodic orbits that are quite close to the network; most likely there are 
unstable periodic orbits as well, but we have not explored this possibility. We 
conjecture that the measure of initial conditions that do not eventually go to 
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Figure 10. Example III: Poincare section, showing five trajectories separated by horizontal lines. 
With xi(0) = 0.01 = h, x 2 {0) = 1, x 3 (0) = y 3 (0) = 10~ 200 , the trajectory starts at the bottom of 
the figure and moves away from the network, but finds a chaotic attractor with < 10 — 67 (below 
the lowest line). This trajectory chaotically switches between ±Y and ±X . Above this, there are 
four examples of stable periodic orbits. The inset enlarges the boxed region near 63 = 7r/2, 
r:j = 10 — 70 . Symbols and colours are as in FigureO 

the network tends to zero as these get closer to the network, so the network 
will be essentially asymptotically stable. The reason for this is that the cusps 
delimiting trajectories that go to ±Y (and all their preimages) are thin. 

4.3 Example III: Trajectories leave the unstable network 

The final example has 5 mm = Sx = 0.99 and <5 max = 5y = 1-10, and shows 
that although the network is unstable, there are nearby periodic and chaotic 
orbits. The Poincare section in Figure [T0l shows five trajectories. Each initially 
has xi(0) = 0.01 = h and X2(0) = 1, and the x-j and ys initial conditions vary 
between trajectories: from bottom to top in the figure, the initial conditions 
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are a? 3 (0) = y 3 (0) = 1(T 200 , 10~ 50 , 10~ 40 , 10" 30 and 1(T 20 , with the resulting 
trajectories being separated by horizontal lines in the figure. Trajectories start- 
ing with the four largest initial conditions find stable periodic orbits; there are 
most likely additional orbits that we have not found. In contrast, trajectories 
starting very close to the network move away and find a chaotic attractor at 
around r 3 ~ 10~ 70 . Within this attractor, the trajectory mostly visits ±X 
and moves gradually away from the network, since 5x < 1. As it does so, 
the chance of visiting ±Y increases; when this happens (with 5y > 1), the 
trajectory jumps closer to the network. The same chaotic attractor was found 
with several initial conditions in the range 

10 -200 < 3.3(0) = 2/3(0) < 1(T 60 . 
This example of sustained chaotic switching between ±X and ±V illustrates 
the ideas discussed in section 13.21 

5 Conclusions 

In this paper we have illustrated a simple mechanism that produces switching 
between different structurally stable heteroclinic cycles in a heteroclinic net- 
work, namely the presence of complex eigenvalues in the linearisation about 
one of the equilibria common to all cycles in the network. This is done in 
the context of an example in M 4 with Z| symmetry. Switching arising from 
the presence of complex eigenvalues has been seen in other examples [1, 2] 
but in those cases the cycles are structurally stable because of transversal in- 
tersections of some manifolds rather than purely because of the presence of 
symmetry. 

By the construction and analysis of maps that model the dynamics near 
cycles in our network, we have found a simple condition under which the 
heteroclinic network is asymptotically stable. The construction of the maps 
used standard techniques that were modified to allow us to keep track of the 
continuum of heteroclinic cycles present in our network. A crucial step in the 
analysis of the network was recognition that the network could efficiently be 
treated as a collection of connections between equilibria (A and B) and an 
invariant circle (C) rather than a collection of cycles each of which connected 
a set of equilibria. Earlier attempts to treat the heteroclinic connections going 
to each of the equilibria on C (i.e., to iX, ±Y,±P and ±Q) separately proved 
to be intractable and were ultimately unfruitful. 
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We found that the network is asymptotically stable if <5 mm > 1 . The quantity 
S mm was defined in section 13.11 and is the product of the ratios of the (real 
part of the) contracting and expanding eigenvalues at selected equilibria in 
the network, namely at A, B, and at the equilibrium on C at which the 
ratio of contracting to expanding eigenvalues is minimised. At B there are 
two expanding eigenvalues, and 5 mm is defined using the larger of these two 
eigenvalues. Effectively, 5 mm is the minimum ratio of contracting to expanding 
eigenvalues that could be encountered by a trajectory on one circuit through 
the network, starting and finishing at one of the common equilibria A or B. 

Another important quantity for network stability is <5 max . As defined in 
section 13. lj, J max is effectively the maximum ratio of contracting to expanding 
eigenvalues that could be encountered by a trajectory on one circuit through 
the network. We showed that if (5 max < 1, then a trajectory started close to 
the network (but not on the stable manifold of any of the equilibria) will be 
further away from the network after one circuit of the network, regardless of 
its itinerary, and thus the network is unstable. 

These results on network stability are a natural generalisation of established 
stability results for heteroclinic cycles, where it has been shown that asymp- 
totic stability of a cycle is often determined (in part, at least) by the ratio of 
contracting to expanding eigenvalues along the cycle. 

We have shown that switching is ubiquitous near our network. In particu- 
lar, we showed that close enough to the network, there are trajectories that, 
over the course of two cycles around the network, visit any combination of 
the equilibrium points within C in any order. (A similar result holds for the 
examples of [1,2].) This occurs regardless of whether or not the network is 
asymptotically stable. In the case that the network is asymptotically stable, 
we showed that most trajectories repeatedly visit both X and —X as they 
approach the network, while, on the assumption that the complex eigenval- 
ues at A mix trajectories effectively, visits to ±Y become rare. Additive noise 
could clearly have an important effect on the switching behaviour; we have 
not explored this issue. 

Our results about repeated switching in the network and our categorization 
of network stability in terms of <5 mm and S ma,x has allowed us to identify an 
interesting case, intermediate between asymptotic stability and complete in- 
stability of the network. Specifically, if 5 mm < 1 < S max then whether or not 
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an individual trajectory approaches the network or diverges from it depends 
on the detailed itinerary of that trajectory. We found two main cases. First, if 
the ratio of contracting to expanding eigenvalues encountered by a trajectory 
making a circuit near the dominant cycle in the network (i.e., a cycle involving 
a visit to either X or —X) is greater than one (5x > 1 ; using the notation 
of section 0]), then we conjecture that almost all trajectories will eventually 
converge to the network, even though the network is asymptotically unstable. 
If on the other hand, the dominant cycle has the appropriate ratio of eigenval- 
ues less than one (5x < 1) with one of the other cycles having a ratio greater 
than one (5y > 1, in the notation of section [4]) then there might be a deli- 
cate balance between the repelling properties of the dominant cycle and the 
attracting properties of the other cycle; most trajectories would be repelled 
from the unstable network but may be attracted to a chaotic or periodic at- 
tractor some small distance away from the network. In this case, the network 
structure may still be observed in the long term dynamics of the system even 
though the network is unstable. 

An example of this phenomenon was shown in section 0] where results from 
numerical integration of a particular system of differential equations were pre- 
sented. We have not attempted to quantify the balance that occurs in the 
example of a chaotic attractor in section HI but defer this to a later paper. We 
note that heteroclinic networks with delicate stability properties have been 
studied before (e.g., in [21]); the point of difference here is that the switching 
mechanism operating in our network ensures that most trajectories will visit 
most parts of a neighbourhood of the network. The transition from 5 mm > 1 
to <5 mm < 1 is an example of a resonance of the heteroclinic network, and it is 
clear from this example that the network structure will make analysis of the 
resonance quite involved. We defer this analysis to a future paper. 
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